/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "distribution.h"
#include <math.h>

static const float halflogtwopi = 0.91893853320467274178;

static inline float
square (float x)
{
  return x * x;
}

static float 
log_normal_pdf (float                   x,
                const Distribution*     ptr)
{
  const NormalDistribution* ndist = (const NormalDistribution *) ptr;

  return - 0.5 * square (x - ndist->mean) / ndist->sq_stddev 
         - halflogtwopi
         - ndist->log_stddev;
}

static float 
dlog_normal_pdf (float                  x,
                 const Distribution*    ptr)
{
  const NormalDistribution* ndist = (const NormalDistribution *) ptr;

  return (ndist->mean - x) / ndist->sq_stddev;
}

static float 
ddlog_normal_pdf (float                  x,
                  const Distribution*    ptr)
{
  const NormalDistribution* ndist = (const NormalDistribution *) ptr;

  (void) x;

  return -1.0 / ndist->sq_stddev;
}

NormalDistribution
normal_distribution (float mean,
                     float stddev)
{
  NormalDistribution rv = { { log_normal_pdf, dlog_normal_pdf, ddlog_normal_pdf },
                            mean,
                            square (stddev),
                            log (fabs (stddev)) };

  return rv;
}

static float
parameterized_log_normal_pdf (float                             x,
                              float                             y,
                              const ParameterizedDistribution*  ptr)
{
  const ParameterizedNormalDistribution* ndist = 
    (const ParameterizedNormalDistribution*) ptr;

  return - 0.5 * square (x - y) / ndist->sq_stddev 
         - halflogtwopi
         - ndist->log_stddev;
}

static float
parameterized_d0log_normal_pdf (float                            x,
                                float                            y,
                                const ParameterizedDistribution* ptr)
{
  const ParameterizedNormalDistribution* ndist = 
    (const ParameterizedNormalDistribution*) ptr;

  return (y - x) / ndist->sq_stddev;
}

static float
parameterized_d1log_normal_pdf (float                            x,
                                float                            y,
                                const ParameterizedDistribution* ptr)
{
  const ParameterizedNormalDistribution* ndist = 
    (const ParameterizedNormalDistribution*) ptr;

  return (x - y) / ndist->sq_stddev;
}

ParameterizedNormalDistribution
parameterized_normal_distribution (float stddev)
{
  ParameterizedNormalDistribution rv =
    { { parameterized_log_normal_pdf,
        parameterized_d0log_normal_pdf,
        parameterized_d1log_normal_pdf },
      square (stddev),
      log (fabs (stddev))
    };

  return rv;
}
